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This paper describes a simple finite element model for large-strain elastostatics. The realization of the model 
in a small-scale computer-code is described. The purpose of the model is to produce test problems for research 
on the application of penalty techniques in nonlinear elasticity. For this reason the code must balance the 
requirements of reasonable flexibility with those of computational economy. The current code employs multilin- 
ear isoparametric elements. The model is capable of generalization to a variety of element types. The solution 
method employed is that of incremental loading combined with the Newton-Raphson method. Symmetric, 
banded systems of equations are produced which are solved in-core. Two- and three-dimensional symmetric 
bodies which are isoparametric images of a reference "brick" may be modeled. An example comparing two- and 
three-dimensional models of a "dogbone" — shaped A.S.T.M. rubber tensile-test specimen is presented. The 
results shed some light on the nature of stress-concentrations which occur in specimens of this geometry. 
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1. Introduction 

The finite element solution of the equations of elastostatic equilibrium for incompressible bodies poses 
special problems not encountered in the compressible case. These problems devolve from the fact that the 
constraint of material incompressibility turns the problem into a constrained minimization problem. There 
is a long history of attempts to find computationally efficient methods to impose the constraint within the 
context of the finite element method [1-5]. ' One of the most promising approaches has proved to be the 
"penalty method" [4-14]. This technique is one in which the discrete equations of the finite element 
method are modified to be able to incorporate a very large bulk modulus. Methods were devised [3] by which 
this could be done without incurring the disastrous loss of accuracy which accompanied early attempts to 
allow large bulk moduli in the finite element equations [1]. The state of the practice of penalty methods is 
now refined enough to allow the bulk modulus to be so large that practically any level of compressiblity can 
be achieved in this fashion. Many problems in incompressible elasticity and fluid dynamics are being suc- 
cessfully attacked using penalties [4-14]. 

The theory of penalty methods has not, however, caught up with practice. The stability and accuracy of 
finite element penalty formulations depend crucially on the choice of element types and numerical integra- 
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tion formulas [3-8], Mathematical theories which can give predicted convergence rates and stability criteria 
exist for problems which are linear or those for which at least the constraint equations are linear [13,14]. 
Most of the successful applications of penalty methods are in such problems [6,8,11,12,14]. When the con- 
straint is nonlinear — such as the constraint of material incompressibility under finite strain — the situation is 
quite different. Therefore problems of large-strain, incompressible elasticity are currently under intense 
scrutiny both in theory and in numerical experiments [9,10,15-19]. 

The purpose of this paper is to describe a relatively simple finite element model designed for the investi- 
gation of penalty methods in finite elasticity. The idea in development of the model is to balance the need 
for sufficient generality with the need for computational economy. Problems with moderately distorted 
geometry may be solved so that the interaction between isoparametric transformations and various numer- 
ical integration schemes may be observed. On the other hand, since it is desired to run numerous test cases 
over and over again varying key parameters, the model incorporates only the simplest of isoparametric mesh 
descriptions. The model is realized in a small-scale code in which all equations can be solved in-core. 

The most important feature of the current model is that it has parallel two- and three-dimensional ver- 
sions. This is designed to confront one of the major obstacles to numerical experimentation in finite elastic- 
ity: the unavailability of exact solutions for all but the most simplified problems. The current model allows 
computation of solutions for thin bodies using a fully three-dimensional formulation. This can be compared 
with plane-stress computations for the same body. In the plane-stress case, the incompressibility is imposed 
exactly in the model. In addition, the plane-stress equations assume a linear displacement variation through 
the thickness of the body. This can be matched exactly by the 3-D model using but one element through the 
thickness. Therefore the degree of success of the penalty method in imposing incompressibility can be 
evaluated by comparsion of the 2- and 3-D results. In this paper, we will deal with a simple, but interesting 
numerical test problem which illustrates such a comparison. 



2. Finite element formulation 

The kind of material properties we shall incorporate in our model are those of rubber-like, hyperelastic, 
isotropic materials. We adopt a total Lagrangian description similar to [19]. The static equilibria of such 
idealized materials are solutions to variational boundary-value problems similar to those described in [20]. 
We have found that we can derive the finite element equations of equilibrium very compactly if we base our 
derivations on differentiation with respect to displacement-gradients, as opposed to Green's strains [19,21]. 
For this reason we are naturally led to the use of the Lagrange stress tensor in our constitutive equation. 

2.1. The variational formulation 

Let u, be the displacement field of a body with reference configuration fl„ in a Cartesian system with coor- 
dinates, x l . If U/j = (diijldxj) is the displacement-gradient tensor, we shall make use of the deformation- 
gradient tensor 



Ju = *v + "u <*> 



and the right Cauchy-Green tensor 



G„ = LJ kJ (2) 

The constitutive equation for the materials under consideration will be of the form 

T = ML = ML (3) 

djji duj,i 

where t is the Lagrange stress tensor and U an energy per unit volume in the undeformed state. 

The energy density, U, is assumed to be expressed as a function of the three principal invariants of Cry 1* 
we denote by ^ y the matrix of cofactors, 
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y s CoUo) (4) 

the three principal invariants of G tJ are given by: 

I = Trace [G y ] = J k J ki 

11 = Trace [Co(G$ = M u (5) 

111= \Gu\ = |/ y | 2 

where | • | denotes the determinant. 

We consider problems whose solutions are characterized by the minima of the variational principle 

*-(".) = Jo. [WWII) - g i u,]dx 1 dx 2 dx 3 - \ T .J t u t ds (6) 

where g t is the body-force and/ the surface tractions, both referred to the undeformed state. The surface of 
I2 , denoted T„, is divided into a load-surface T°and a surface with prescribed displacements, F°. All loads 
are assumed to be dead loads. 

2.2. Isoparametric finite elements 

We shall consider a particular kind of finite element mesh. It will be composed of only one type of ele- 
ment, which is defined on a reference domain Q R (the unit cube, for example), and is mapped into its loca- 
tion in the global mesh by a transformation based on the element shape-functions defined on £l R . A variety 
of moderately complex bodies may be modeled in this fashion, and with the development of some simple 
notation, the model can be described without much difficulty. 

The process of constructing a finite element assemblage involves piecing together small domains to form 
a large domain. Similarly, the global trial functions are the unions of element trial functions which are sup- 
ported (i.e., are nonzero) only on individual elements. The global trial functions we describe are compatible 
for elasticity problems [22,23]. This means that they are continuous across interelement boundaries, which is 
assured by requiring that the values of the trial function agree on shared boundaries. 

o. The element transformations 

Let the i ,h component of the local coordinates of Sl R be £', and select & as coordinates of the n' h of n R 
nodal points of the prototypical element. Construct n R shape functions, <f> n defined for all values of £'. These 
shape functions are polynomials on Q R , such that 

<M£') = ?<Q K 

(7) 
<M£) = S„ r,s = l, . . . ,n R 

according to the procedure outlined in [22,23] and illustrated by an example in figure 1. 
For each of the 1 < e < N e elements in the global mesh, define a vector-valued map 

T? a t\ r 4> r (8) 

The summation convention is used on the non-tensorial index r = 1, . . . ,n R , but is never used on the index 
e. Equation (8) defines the isoparametric transformation of each element from fl* to its location in the global 
mesh. Note that the global nodal locations are then 

TM) = t% (9) 

Examples of typical transformations are illustrated in figure 2. 

We assume that the t' r are chosen to satisfy certain basic properties: (1) that compatibility for elasticity 
problems is satisfied by the global shape functions, (2) that the totality of the transformations ^'defines an 
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r=4 



{' 



FIGURE 1. A finite element shape function on a reference four-node element. 




(-1.-D 



Figure 2. The transformation by T' from the reference element on 0« to the global mesh on Q„. 

assemblage of elements such that their domain Q with boundary F~° approximates ft well enough so that U a 
— 0„ is negligible for the purposes of the model, (3) that J) is an invertible transformation at every point of 
element e, with inverse S*. The inverse 5* is extended to o by 



Si&) 



-{ 



£', 9 Tj(£') = x>, x> e element e 
any point not in Q*, otherwise 



(10) 



The meshes which we will employ in our model are simple enough that the connectivity of the global mesh 
can be described by a closed form, integer-valued function, A(e,r), called the "nodal alignment function". 
To each element e, 1 < e < N. and r, 1 < r < n R , A assigns the global ordering number of the node with coor- 
dinates t ir : For more complicated meshes, A is well-defined, but may be determined from a table of data, 
rather than an arithmetic expression. In terms o(A(e,r), the global shape- functions have the form 

*' = U (0 r [S;(.)]} 

9/t(,,r)=/ 



(ID 
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In terms of (11), it can be seen that the purpose of the extension of S.'defined in (10) is merely to assure via 
(7) that whenever^ is not in any of the elements with A(e,r) — I, then <t>X#0 is well-defined and equal to zero. 
A typical $»/ and its construction is illustrated in figure 3. The finite element trial space, S\ is then the 
space of all trial functions which satisfy the displacement boundary conditions on T* and have the form 



Ui(x J ) = u't^iix 1 ) 



(12) 



The sum on / is taken over all 1 < / < N global ordering numbers, and u\ is the i ,h component of the nodal- 
value vector at node /. Equations (7) and (11) imply that if x{ are the coordinates of the /'* node in the global 
ordering, 



Ui{x{) = u\ 



(13) 




A (34.4) = 50 
A (35,3) = 50 
A (37.2) = 50 
A (38.1) s 50 



FIGURE 3. Alignment of local and global numberings, showing the role of A(e,r) in piecing 
together the global shape function 4> so from the <j>, on fl*. 



The restriction of u.- to element e can be determined from the vectors of local nodal values in element e, u", 
r = 1, . . . , n R : 



uV = u\, A(e,r) = I 



(14) 



Then by (11) through (14), the restrictions of u t and u u to element e, denoted by, ufand u'j, as functions of 
the reference coordinates, are 



uijf) = [{dT-ruV -|^r Oj] (£*) 
ok 



(15) 



where £* = 5j (*>) and dT* and Q, are defined by 
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dT< = \dn\ (i6) 

Qj = Co{dT'j) 

To avoid confusion, differentiation with respect to the x' will use the comma notation, and differentiation 
with respect to the £' will use dld£'. 

b. The element energy expressions 

For U evaluated at any u, e S h 

|jj o Udx l dx 2 dx 3 = L. 1 0r U e dT'd?d?d? (17) 

where 

U. = 0(1., II., III.) 

(18) 
I. = I [</(£*)], etc. 

The integral is evaluated by a quadrature rule with weights w k and evaluation points rji on £2 R . We define 
elemental weights on e by 

wl = w k dT'(nQ (19) 

The final energy expression is then 

J U dx l dx 2 dx 3 = L, L k wl U e (7/i) (20) 

c. Modifications for penalty methods 

The numerical integration scheme of (20) is adequate for a compressible energy density. For incompressi- 
ble or nearly incompressible energy densities, following [8-10,15,19], we have 

tf(I,II,III) = »UII) + y G 2 (III-l) (21) 

where z is a bulk modulus scale factor; z may also be interpreted as the penalty enforcing III-l to be arbi- 
trarily close to zero for material incompressibility [8-10,15,19]. In the plane-stress case J u is assumed to have 
a form which enforces III = 1 so that £/(I,II,III) = W(\,\Y). This will be fully described in section 3.2. It 
allows the numerical integration scheme of (20) to be applied in the plane-stress case. 

For the 3-D model, two integration formulas are required, the one in (20) applied to the W terms in (21), 
and a "reduced" formula, with weights r,and evaluation points ft' applied to the G 2 terms [8-10,22]. Equa- 
tion 19 has analogue: 

r' ( = r ( dT'(M (22) 

and (20) is replaced by: 

J Oo U dx'dx'dx 3 = E. [E k uf k WM) + f£, >1G 2 (ft)] (23) 

where W, and G, are defined as in (18). 
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2.3. Derivatives of the energy-density 

The load terms in (6) from body and surface forces may be treated exactly as in linear finite element 
analysis [22,23]. Combining this with (20) or (23) gives a numerically integrated variational principle to be 
minimized over S h . We will denote this principle by the same 7r(u,) and understand that throughout the 
following discussion the minimum is to be taken over S\ The equations of equilibrium of the finite element 
model are then 



(24) 






i = 1, 2 and / = 1, 2, . . . , N 
in the plane stress case, and 






+ -f- E f-^r[«](W] +/? = (25) 

» = I, 2, 3 and / = 1, 2, . . . , N 



For the 3-D case. 



f\= £ //" (26) 

9^(.,r) = / 

where f" is the i"* component of the distributed nodal force from body and surface forces acting on node r of 
element e [22,23]. Thus /' is just the usual finite element load vector, and (24-26) express the fact that the 
assembly of the equilibrium equations of the nonlinear model obeys the same rule as assembly of load vec- 
tors in linear analysis [22,23]. 

The equations of equilibrium, Ff = 0, will be solved by a combination of incremental loading with the 
Newton-Raphson method [22]. More details of these procedures will be described in section 4. Both incre- 
mental loading and the Newton-Raphson method require calculation of the "tangent stiffness matrix" [22], 

K, = [^] (27) 

Since K r is the matrix of second derivatives with respect to u' of tt, which is a continuously differentiate 
function of the nodal- values, K r is a symmetric "hessian" matrix. 

"atf-isr - £> LkWk duvduy (,J (28) 



in the plane-stress case and 



W - dFj _ T r r . dW. , 



In other words, (27-29) imply that K r assembles by the same rule by which stiffness matrices assemble in 
linear analysis [22,23]. Beyond that, K r will have the same banded structure which a stiffness matrix in 



duVdu'j 
in the 3-D case. 
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linear analysis would have, using the same mesh [22,23]. This follows because the band-structure is deter- 
mined by 

| I-J\ (30) 



K r will have (2oj + l)d bands, where d = 3 for the 3-D model and d =2 for plane-stress. Clearly w is deter- 
mined by the form of A and not the form of the energy expression. Boundary conditions are also imposed on 
V u in the same manner as in linear analysis [22,23]. 

a. Separation of tensorial factors 

In what follows, let U, represent generically either a compressible element energy density, (7(I.,II„III.), or 
one of the terms of an incompressible element energy density, W. or G t . From the previous section we may 
conclude that 

dU, _ dU. du' Kt _ dU e t ^ T . Y i d<fr r r , 

(31) 
d 2 U. _ d 2 U. ( . Tty2 d<f>, dct>, 

faVduV ~ Bu' u du' k , t (dI > d£ m d? LmjLnt 

We point out again that there is no summation on e; in fact, since in what follows we will always be con- 
cerned with elemental expressions, we will drop the index e. At this point the strain-energy derivatives are 
seen to be the sum of terms which are separated into factors, the second of which, 

(W -fgr C mJ or 077* jfr -ffr C^ (32) 

is determined from the elemental isoparametrics and local shape functions, and is no different in character 
from the corresponding factors in linear isoparametric problems [22,23]. What distinguishes the finite 
elasticity case is the complexity of the first factors. It is these derivatives of the energy density with respect 
to displacement-gradients to which we now turn our attention. Actually, we may further restrict our atten- 
tion to the derivatives of strain-invariants, since 

U, = ML = MLM- + MLML + dU dill 

du t j 31 du t j 311 dutj dill du tJ 

(33) 

U = d 2 U = dHJ__dl dl_ + dU d 2 l 

duijduk.i dl 2 duu du kil dl du,ijdu k<( 

+ similar terms for II and III 

au a 2 T7 an 

,. > - T2 > „. • • . , etc., depend on the form of the energy density for a specific material and are not 
usually difficult to evaluate. So the problem is to compute the derivatives of the invariants. 

b. Derivatives of the strain-invariants 

Derivatives of the strain invariants, I y , I w , II y , . . . , etc., are defined in a fashion similar to U u and U m . 
The expression for I v follows immediately from equation (5), as does III y , by remembering that the 
derivative of a determinant with respect to the matrix itself is the matrix of cofactors. Next we observe that 
the derivative of the matrix of cofactors with respect to the matrix itself can be written 

"a"r = £ikmCjtnJmn (34) 



lk t 
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where e tJk is the Cartesian permutation symbol. Using the chain-rule and collecting terms by application of 
symmetries gives an expression for Ily. So we have 

Iy = 2/ y 

Hy = 2eikm€j(nJmnYkl /gin 

m v = 2\j v \* u 

Calculating I w is trivial. Using the chain-rule and the identity e Up € kt p = 8 lk 8 J( — 5/, 5 V leads readily to an 
expression for Hy*,. To find Illy*,, use can be made of the inverse of /y, (/y)" 1 = $/,/ |/*f|. Then Illy 
becomes 

Illu = 2\J u \Vj,r (36) 

(note the transposition). Making use of the identity 

-^P- = - (/„)" W = Mvl |/« | 2 (37) 

leads to a simple expression for Illy*,. Accordingly we have 

Iy*f = 25ikOj t 

\\ ukt = 2e ikm e Jln \fr mn + 28 ik 6 Jt J mH J mn (38) 

— 28 ik J m( J mJ — 28 Jt J im J km + 2JtjJ kl 

Illy*, = 4^y0*, - 2\l/,f\l/ kJ 

Although the terms with the e iJk may look forbidding, they may be computed using no floating point opera- 
tions other than "change sign," and the same calculations can be used for Hy and Ily*,. Alternative expres- 
sions (of comparable complexity) for Hy and Hy*, can be derived without e y * using the identity: II = 

(l/2XF-GyGy). 



3. Two- and three-dimensional problems 

In this section we describe the assumptions which are made in reducing a three-dimensional problem in 
finite elasticity to a two-dimensional problem. We describe how this reduction is incorporatd in the 
framework of the results of section 2. We also will describe the similarity and differences between the ele- 
ment shape-functions and the isoparametric transformations in the two- and three-dimensional cases. It will 
be easier to describe the implementation of the isoparametric transformations in two-dimensions. The three- 
dimensional implementation can be viewed as a straight-forward generalization of the 2-D case. Finally, as 
will be explained in more detail in section 5, linear elastic problems serve an important purpose in computer 
code verification. They are also interesting in their own right. We therefore conclude this section with a 
discussion of two- and three-dimensional linear elasticity problems. The linear two-dimensional model incor- 
porates compressible, nearly incompressible and incompressible materials in states of plane stress or plane 
strain. Penalties are used to enforce plane-strain incompressibility. 

3.1. Three-dimensional problems 

The results of section 2 were derived in full three-dimensional form. They may be applied directly to the 
coding of the element matrix, element load-vector and assembly subroutines [22,23]. Details of how the mesh 
is described and how the isoparametric transformations are generated will be described later in this section. 
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But the two-dimensional, plane-stress model requires additional development. It may be thought of as a 
three-dimensional model in which specific simplifying assumptions are imposed on the stress and strain ten- 
sors. To see how this is accomplished, it is useful to look first at those tensors in the fully three-dimensional 
case and observe ah example of a particular material model. 

a. The general Incompressible formulation 

In the 3-D case J u , Gy, T y and all the tensors of the previous section are of dimension three in each rank. 
The Lagrange stress tensor can be written in terms of (35) 

For incompressible, 3-D problems, U is assumed to have the form (21), so (39) becomes 

T " = "^P lji + -flf H " + zG ( Ul ~V Jfr HI* (40) 

Recall that z is a large penalty parameter forcing G(III-l) to be small in the finite element solution. 
Therefore the pressure is given by 

h=zG{\\\-\)-ML (41) 

so that (40) gives an incompressible constitutive equation in the z— 00 limit. 

For stress calculations from the finite element model, it is more physically meaningful to compute the 
Euler stress tensor, referred to the strained body: 

fry = /*,T w /|/ mn | (42) 

Equation (35) and the fact that (Jj)' 1 = ^y/|/y| with (41-42) imply 

olj = [-^P M* + -|^ ] k \\ Jk + 2MI5.J/ |/ mn | (43) 

Then denote by B tJ the left Cauchy-Green tensor and by C v its inverse 

Bij = JikJjk 

(44) 
Cij = MJ\J mn \* 

We can rewrite II as 

II = C u \G mn \ (45) 

from which it easily follows using (37) that 

Jk,llj k = 2lll[C kk 8 iJ - C tJ ] (46) 

which leads to 

a; = 2Vm[-flB t ,--|^ c „ + /;5j (47) 

in which 

88 



h = nih + IK- Clt (48) 

oil 

Finally it is desired to have a stress-free undeformed state, but for J u - 5 U (47) gives 

*; = AA, = (4f +2-f{fR (49) 

so letting H = h — h 

a, = 2Vm [&- B tJ - -§jf C„ + H5,] (50) 

All constituents of (50) are easy to compute given J v . Symmetry implies only six of the components of Oy 
need be calculated. 

b. Application to Mooney-Rivlin materials 

A simple constitutive equation of the form (40) can be obtained by taking 

dl " dli Cl 

(51) 
G(x) = x 

This is a penalty function version of the Mooney-Rivlin constitutive equation [9,10,19, and 21]; Ci and c 2 are 
material constants. This equation is intended to model exactly incompressible materials, and in fact, it has 
some continuum mechanical inadequacies if significant compressive deformation is allowed [10]. Therefore 
it should only be employed with large values of z. The exactly incompressible, plane-stress version of the 
constitutive model has been found to give reasonable agreement with experimental results in some problems 
[21]. How well the 3-D penalty version succeeds will be discussed below. 
The symmetric stress tensor for Mooney-Rivlin material is given by 

o tj = 2yfm [ Cl B u - c 2 C y + Hb tJ ] 

(52) 
H = 2(111-1)111 + c 2 C i{ - c, - 2c 2 

It should be noted that this differs slightly from the /fused elsewhere [9,19,21], because c 2 Ciib K i j is considered 
a pressure term. This means that H defined in (52) has the value c 2 -Ci in the undeformed state. [In [19], 
since the CiCubij term is included with the deviatoric terms, the pressure in the undeformed state is 
— Ci— 2c 2 .] A quantity of more interest than H, then, is the pressure above the undeformed state: 

H = 2 (III - 1)111 + c 2 (C„-3) (53) 



3.2. Finite plane stress 

a. The general incompressible formulation 

The assumptions in the plane-stress model are that 

1. The body is very thin, with uniform thickness, 2t, in the r'-coordinate. How thin the body must be 
will be dealt with further below. 

2. The applied forces are in the x l -x 2 plane; the upper and lower surfaces in the x x -x 2 plane are free 
surfaces. 

3. a 3l — o 32 = o 33 = throughout the thickness, i.e. for —t :< x 3 < t. 
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4. The deformation is essentially biaxial, i.e. 



u 2 {x\x 2 ,x y ) = u 2 (x 1 ,x 2 ) 
^(x'.x 2 ,* 3 ) = [9(x',x 2 )-l]x 3 



(54) 



for an unknown function ©(xSx 2 ). 
5. The body is exactly incompressible. 



This model is discussed in more detail in [21]. An example of such a deformation is shown schematically in 
figure 4. 




u 



u_(x , .x 2 ,x 3 ) = Ds , (x , ,x 2 )-l]x 3 



FicURE 4. Plane-stress assumptions for the deformation through the thickness of a thin body. 
These assumptions imply that the deformation-gradients have the form 



Ju = 



"1+",,, 


"1,2 


0' 


"2,1 


l+"2,2 





^ * 3 e,, 


x 3 e, 2 


e 



(55) 



Assumption 1 is further specified by the requirement that the body be thin enough for tQ A and tQ^ to be 
neglected. Note that this is not an absolute requirement on the thickness of the body, since if 0,i and 0, 2 are 
small enough, it can compensate for a larger /. On the other hand, very rapid spatial variations of Q{x l ,x 2 ) 
could invalidate the assumption, even for small t. 



Let K ag be defined for a,/3 = 1,2 



Arf — 



l+"l,l "lj 

L "2,1 1+"2J 



(56) 



We will use Greek subscripts to indicate tensors which are of dimension 2 in each rank. Imposing the thin- 
ness assumption leads to 

Ju = [K* © ©J (57) 

where " e " means "matrix direct sum" of the 2x2 matrix, K a/3 , with the lxl matrix, 0. That is 

K a g 



Ju = 











(58) 
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Assumption 5 is enforced by taking = A~\x l ,x 2 ) where 

A(x\x 2 ) = \K a0 \ 



(59) 



Now we can rewrite the invariants of Gy and their derivatives, using (35-38) and the special form (58). In (58) 
we have mixed tensor and matrix notation in an obvious way. In what follows, we wish to construct tensors of 
rank 4 from tensors of rank 2 given in matrix notation. We adopt the convention that a matrix with 
subscripts appended on the lower right defines a tensor of rank 2 with components given by the matrix 
entries. The first subscript is the row index and the second the column index, as usual. For example, in 
terms of (58), ip u of cofactors of J v is 



Similarly 



V = [[A- l CoK a0 ] © A}„. (60) 

I = K a0 K a0 + A' 2 

II = A- 2 CoK a0 CoK a0 + A 2 

I„ = 2[K a0 © A" 1 ],, 

II y = 2[[A- 2 K ag + ACoK a0 ] © [A- l CoK ag CoK a0 ]}u 

Illy = 2{[A- 1 CoK a0 ] © A} l7 (61) 

lijkt = 25 ik 8jt 

U ukt = 2[e ikm eji n ipmn + 5 ik djtl 

- 8 ik [K ya K y0 <BA- 2 ] Jt - dAK^QA- 2 ]* 

+ JuJki]} 

Equation (58) has made III = 1 and the energy density is not a function of III, so IIIy*f is not needed. Illy 
would only be needed in the event one wished to calculate r y . By enforcing plane-stress assumption 3, the 
pressure can be determined directly. It should be observed that (58) and (59) already imply that a 3 i = a 32 = 
0, since (58) implies 



B iJ = [K ay K 0y eA- 2 ] IJ 
C tJ = [A- 2 CoK ya CoK y0 eA 2 ] tJ 
Substituting (62) in (50), setting a 3 3 = and solving for H implies 

H= -lf. A - 2 +-^-A 2 

di an 

The reduced constitutive law may be now written 

, U = 2[&-B u -*3LC V + Hii 

H = -M-A- 2 + dW A2 



(62) 



(63) 



(64) 



di 



dll 
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Assumption 3 is enforced by choosing plane finite elements of the type illustrated in figures 1-3 and by 
replacing every occurrence of u„, 3 or u 3>a by zero (or = 1,2) and u 3J by in the energy expressions of section 
2.2. Any explicit occurrence of* 3 is replaced by x 3 = 0, since we assume negligible variation of all quantities 
through the thickness. The ac 3 -integral in the multiple integrals of the energy expressions integrates to a fac- 
tor of 2t in the energy expressions. This may be ignored since this will not change the minimizing displace- 
ment field. The result of the reduction is that there is no dependence on x 3 or t in the plane-stress model. 

A reduced displacement field is used which has only two components, u a ; u 3 is a function of u x and u 2 . 
Therefore in section 2.3, all derivatives with respect to nodal values of u 3 vanish. But it is important to note 

that in (31), terms of the form „ ' do not vanish. They have a new form, introduced by the plane-stress 

OUa 

assumption of the form of u 3 (x l ,x 2 ). For notational convenience, define 

(t>r 



-orr-ffe-c*, 



d? 



(65) 



This is the r' h element shape function's derivative with respect to Cartesian spatial variable x m , and for 
plane stress m = 1 or 2. Let L = I, II or III and u r a a = 1 or 2 be the <x' h component of the r" 1 local nodal 
value in an element 






dL 
du r a 

d 2 L 

duadu'g duijdu k ,i du r a du% 



d 2 L 



du' c 
diijj du k j 



+ L U 



d 2 u u 

dUadu'g 



(66) 



- T toy dUk ' . t d 2 "ij 



dua dui 



duLdue 



Since trial functions are linear in the nodal-value coefficients u r a , if i andy =£ 3 we have 

du 'J _ * . d 2 u tJ _ 
o ia (p rJ , ^ _*\ = 



bu 



du'du'n 



(67) 



u 3 j is not linear in u r a and is a function of the in-plane displacement gradients. Therefore contributions to 
the chain-rule sums in (66) from i,j,k, and t = 1 or 2 can be computed directly from (61), and (67) which hold 
for these values of the subscripts i,j,k, and I The contributions to (66) from i = j = 3 and/or it = f = 3 may 
be found by using L u and L iJkl determined by (61) and the following facts: 



du X3 _ dA~ 
du r a ~ du r a 

3A 



= - A~ 



BA 
du r a 



Q u r a - ( l) a * l {<f>r a J 3 - a 3 -« ~ <f>r 3-c.L-a 

(no sum on a) 



J 



dA 



du% 



dUadll'g QUa 

•^7 = (l-U(-i)-(0 ra ^-0, 
(no sum on <x,p) 

b. Application to Mooney-Rivlln materials 



d 2 A 

dUadll'f} 
S0*») 



When the constitutive equation with ■—- and ^-~ a\ 



(68) 



(64) becomes 



dl 



otj- given as in (51) is employed in the plane-stress case 



°u = 2[c 1 B iJ - c 2 C tJ + Hb tJ ] 
H = -c, A" 2 + c 2 A 2 
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(69) 



H, the pressure above the undeformed state, given by 

H = Cl (l-A- 2 ) + C2 (A 2 -1) 
3.3. Shape functions and isoparametrics 



(70) 



Here we discuss the specific choice of shape functions, <j) r , and the isoparametric transformations, 7), 
which are described generally in section 2.2. It will be easier to describe the two-dimensional case first and 
think of the 3-D case as the appropriate generalization. 

a. Two dimensions 

The reference element, 12*, is the square in the £* — £ 2 plane with — 1<£'<1. On this domain, 



<M£\£ 2 ) = 



44JI 2 



(71) 



where (& #) = ( - 1, - 1), (& £§) = (1, - 1), (& ffl = ( - 1,1) and (& ® = (1,1). These element shape func- 
tions look exactly like the one illustrated in figure 1. 

To construct the isoparametric transformation, T', all that is needed is to determine the nodal locations, 
t' r , of (8). Then the isoparametric transformation of 12* is determined throughout element e via (8). The 
nodal locations f ir give the £ = 1,2 coordinates in the global mesh of the nodes r = 1,2,3,4 of element e. 
From this information, the rest of the computations implementing the ideas of section 2 can be carried out 
once the nodal alignment function, A{e, i), is defined. To do this, we chose to restrict our attention to curved 
bodies which are described as a class of differentiable, invertible tensor transformations of the reference 
"brick" of elements illustrated in figure 5. It is clear that the nodal alignment of such a body is the same as 
the nodal alignment of the reference brick. Therefore the global element and nodal numbering schemes and 
A(e,i) may be determined from the reference brick. We have chosen the numbering scheme illustrated in 
figure 5. In terms of these numberings, A(e,i) can be expressed as a simple, closed form, integer expression. 
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Figure 5. A two-dimensional reference brick before and after curving. Numbers indicate beginning and 
ending of global numbering scheme. 
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As a further simplification, we consider only bodies which have symmetry lines on the lower boundary and 
rightmost boundary as illustrated in figure 5. The transformation of the reference brick is generated by 
defining a cubic spline [24] which interpolates to given data defining the location of the upper boundary of 
the curved body (the leftmost boundary is not curved). This spline is used to define a transformation which 
varies linearly in the interior of the reference brick in such a way that the reference brick and curved body 
have the same axes of symmetry. This means that the coordinates, f ir , of the nodal locations in the curved 
body can be obtained by multiplying the appropriate coordinate of the corresponding nodes in the reference 
brick by a factor which depends on the given cubic spline and the distance of that node from the symmetry 
line in the reference brick. 

b. Three dimensions 



The <j> r are trilinear polynomials on Q R , which is the cube such that -1<£'<1. 

(* , +4;)<! a +«?)G 3 +E) 



<Mr,£ 2 ,£ 3 ) = 



si, 1 !?! 3 



(72) 



where £' are the coordinates of the corners of $l s . 

The f ir are determined analogously to the way they are determined in 2-D. In 3-D, the back, rightmostand 
bottom surfaces of the reference brick are symmetry surfaces, as illustrated in figure 6. The leftmost surface 
is not curved. The remaining two surfaces are defined by interpolating splines and define transformations 
which vary linearly in the interior in such a way that the symmetry planes of reference and curved domain 
coincide. In the current code a further simplification is obtained by requiring that the surface splines be 
tensor products of one-dimensional splines. This restricts the type of surface which can be described but 
allows the use of a simple, one-dimensional spline code. 



/~7^ 
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Figure 6. The three dimensional generalization of the curving pro- 
cedure illustrated in figure 5. 
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3.4. Linear elasticity 

The current code can compute finite element solutions to standard 2- and 3-D linear isotropic elasticity 
problems [22]. The same shape-functions and mesh generation are used as for the nonlinear problems. The 
linear solutions can be compared to small-strain nonlinear solutions. This allows determination of the load 
and/or strain levels at which the linear and nonlinear elastic models depart significantly in their predic- 
tions. It is also interesting to observe the nature of the differences in the linear and nonlinear solutions. A 
simple analysis shows that the Mooney constitutive laws (52) and (69) have shear modulus 

M=2(c, + c 2 ) (73) 

in the zero-strain limit. The incompressibility condition determines the other constitutive parameter of 
linear, isotropic elasticity which give a correspondence between the material properties of the linear and 
nonlinear models in the small-strain limit 

For comparison with nonlinear plane-stress solutions, it is natural to look at linear plane-stress solutions. 
As with the nonlinear constitutive laws, the linear plane-stress model enforces incompressibility exactly by 
assumption, without penalties. However, by a well-known re-interpretation of the material properties, every 
plane-stress solution is also a plane-strain solution [25]. Incompressible plane-strain requires a penalty to 
enforce zero volume change. Therefore the current code can be used to evaluate the behavior of two- 
dimensional elements with penalties, at least in linear elasticity. 

a. Three-dimensional linear problems 

Three-dimensional linear elasticity with a penalty to enforce incompressibility is treated exactly as 
described in [5]. In terms of the infinitesimal strain tensor e u = l/2(u iJ+ u,,,), the constitutive law is 

o y = 2 M {ey + (zeM (74) 

\x is the shear modulus and z is the penalty parameter 

1 —2v 
v is Poisson's ratio, which has the value 1/2 for exact incompressibility. By taking 

h = -J- e„ (76) 

equation (74) becomes 

a iy = 2/i [e u + Ph8 u ) (77) 

which is the usual constitutive law for linear, isotropic elasticity, which is valid for incompressible material if 
v = 1/2 [1]. 



b. Linear plane-stress and plane-strain 

The linear plane-strain constitutive law is used to model the cross-section of a body which is very long in 
the x 3 -direction [25]. The problem is two-dimensional since it is assumed that M.Cr 1 ,* 2 ,* 3 ) = Ui(x l ,x 2 ) and 
^33 = 0. Substituting this into (74-77) gives a 2-D constitutive law for which 

h = -*-(<?„ + e 22 ) (78) 
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Then as ,-1/2, z is a penalty parameter enforcing e u + e 22 = in the incompressible limit. The energy 
expressions and the finite element model are very similar to 3-D linear elasticity, but e 33 - 0, and there is no 
x 3 -dependence. If we take a,j3 = 1,2 the plane-strain constitutive law is 

a ag = 2fi{e a0 + * 5 af ,e yy ] (79) 



For any value of v ^ V2, if we re-interpret a solution obtained using (79) by defining 

(80) 



" = Tf7 



e 33 = - 1 _ 0?ii + e 22 ) 

it is easy to see that we have a linear plane-stress solution satisfying o 33 = 0, with apparent Poisson's ratio n, 
and shear modulus /i [25]. In particular, if we want the plane-stress solution to be incompressible, we choose 
77 = V2. This corresponds to v = Vz in (79), or z = 1. The current code then can be used to obtain an incom- 
pressible plane-stress solution by setting v = Vz, and an incompressible plane-strain solution by the penalty 
approach, taking v close to V2, but not exactly equal to V2. Since they involve quite different geometric 
assumptions, the two kinds of 2-D incompressible solutions usually are quite different. The same kind of 
comparisons between the plane-stress solutions and 3-D solutions are possible as were for nonlinear 
elasticity. 



4. Nonlinear equation-solving procedures 

One of the principles which has guided us in the construction of this finite element model is that it should 
retain as much of the structure of a linear elasticity model as possible. Other than in the assumption of some- 
what simplified loading, this has been done without making any restriction on the degree of nonlinearity in 
the energy expressions. The point is that this can be done because the energy principle (6) which charac- 
terizes solutions in our chosen continuum model admits to the generalized Rayleigh-Ritz procedures 
described above. Because trial function are linear combinations of the basis functions, and because of the 
linearity of the integral in the energy expressions — features which are common to linear and nonlinear 
elastic models — the assembly of elemental expressions and much of the generation of these expressions are 
essentially linear procedures. Much of the current computer code is the same code which would be 
generated to solve linear elasticity problems alone. 

To solve the nonlinear equations (24) or (25), many of the common procedures involve generation of a 
sequence of linearizations using the tangent stiffness matrix (27). A sequence of linear systems is solved with 
K r on the left-hand side. We have seen that K r has the same banded structure in the linear and nonlinear 
case. Furthermore, it is a consequence of variational calculus [26] that at the solution to (24) or (25), K r is 
positive-definite (we assume our problem is well-posed). It is a guiding principle of the nonlinear iteration 
schemes described here to retain positive definiteness of Kr throughout the whole sequence of iterates. First 
of all, this will allow even more common code between the linear and nonlinear procedures: the same stable 
positive-definite, linear-equation solvers can be used. Second — and in our view more important — is that 
departure of the sequence of iterates from the region of positive-definiteness of K r has serious implications 
for the stability of the whole solution process — both in a numerical and analytical sense. Loss of positive 
definiteness of K T means that the numerical inversion procedure is less stable in the sense that small round- 
ing errors are more likely to be magnified [23]. But potentially more serious is the fact that loss of positive- 
definiteness of K r is an indication that the current iterate in the nonlinear scheme is outside the region of 
local convexity surrounding the desired solution [21]. The iteration scheme runs a danger of becoming 
"lost" and failing to converge. 

The numerical results described in the next section show that the possibility of the iteration scheme 
becoming lost in a region of nonconvexity of the energy is sometimes a very real one — in spite of the precau- 
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tions we take. One of the most important conclusions which can be drawn from the current state of research 
into penalty methods is finite elasticity is a negative one: penalties greatly increase the likelihood of depar- 
ture from the convex region. Before we discuss this further we describe the iteration schemes we employ. 

4.1. Incremental loading 

Iteration schemes require an "initial guess" solution at which to begin. The degree of nonlinearity in 
equation (24) or (25) is, to a large extent, dependent on the intensity of the loads as reflected in the 
magnitude of the load terms. A useful strategy in producing an initial guess at the desired load-level is to 
begin from the unloaded state with u, = 0, and increment the load gradually to the desired level. At each 
step of this "incremental loading" process [21,22], an approximate solution is generated which is an initial 
guess for the next load-level. When the final load level is achieved, the approximate solution at that load can 
be used as an initial guess for another iterative method which produces the final, refined solution. The 
strength of this strategy lies in the fact that for small enough load increments, the approximate solution at 
each load stays in a region of convexity of the energy. Furthermore, if the errors which accumulate during 
the loading process tend to cause departure from the region of convexity, the user has three options: (1) 
refine the load increment; (2) use the other iterative method to eliminate the accumulating errors at inter- 
mediate load-levels; or (3) combine (1) and (2). 

Let u represent the nodal-value vector of displacements in (24) or (25). Let the strain energy terms in these 
equations be denoted by the vector P(u) and the load terms by F. If < X < 1 is a scalar, at any stage of the 
incremental loading process we solve the nonlinear equation 

P{u + Au) + (X + A\)F = (81) 

where AXF is the load increment, and Au is the change in the displacement field in moving from solution u 
at load-level XF to solution u + Au at load (X + AX)F. Expanding P(u + Au) in a multivariate Taylor series 
and retaining only first-order terms in Au, along with the definition of K r = K T (u) in (27) leads to 

P(u) + K T (u)Au + XF + AXF = (82) 

If u is the exact solution at the previous load-level or u has been obtained by a previous step of incremental 
loading, P(u) + \F = 0[(AX) 2 ]. This gives an algorithm with 0(AX) accumulated error: 



u <o> =0, Xo = 

X, = X,-, + AX 

(83) 
Au = -AXKf 1 (" (, " ,, )F 

«<••) = «<'-» + Au 



The iterations in (83) can be carried out until A = 1, or interrupted before then so that the current u U) can 
be replaced by a refined u (,) . The algorithms which do the refinement (such as the Newton-Raphson method 
described in the next section) tend to be more costly than one step of incremental loading. Therefore in the 
attempt to remain in the region of convexity, a judicious balance between load-increment, AA, and number 
of intermediate refinements may be called for. In our experience with the Mooney-Rivlin constitutive law, 
convexity can be maintained with very large AA, so that no intermediate refinement is needed, in the 2-D 
plane-stress model. On the other hand, efficient computations in the 3-D case require a careful balance 
between choice of AA and the number and frequency of intermediate refinements. This is due to the sensitiv- 
ity of the incompressibility constraint equations to small volume changes [9,10]. 
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4.2. Newton-Raphson iterations 

Let kF denote the current load-level, for < A < 1. If u is an initial guess at this load-level, obtained either 
by incremental loading or other means, we wish to compute a correction, Au, so that 

P(u + Au) + \F=0 (84) 

Again making a Taylor series expansion for P(u + Au) and neglecting second order terms in Au yields the 
Newton-Raphson method: 



Au = - K? 1 (u ( '- 1 »XP(u (, '- 1) ) + XF\ (85) 

u (.) = u ( '"- l) +Au 

The neglect of second order terms in (85) means that several iterations may be needed. It can be shown that 
in any vector norm ||u (, ' +,, -u (,) || = 0(||u (, ' ) -u (, '~ 1, || 2 ) [27], if u (,_,) is "sufficiently close" to the exact 
solution. This means that convergence can be very rapid. But it also implies that no prior guarantees can be 
made if u u ~ l) is far from the exact solution. The algorithm can get lost in a nonconvex region, even if it 
begins with u in the desired region of convexity. We have found that with the 2-D problems, almost any u 
obtained by incremental loading (even with AX large) is sufficiently close, right from the start of (85), to 
initiate quadratic or near quadratic convergence. Furthermore, even if u is far enough from the solution so 
that the ultimate quadratic convergence rate is not achieved from the start, the iterates of (85) converge 
more slowly until u u_1) is close enough to give quadratic convergence. 

In 3-D problems, we find that care must be taken in choosing AA so that the initial guess, u, is accurate 
enough for (85) not to get lost. If this is done successfully then the ultimate quadratic convergence rate is 
guaranteed. The dominant computational cost in (85) is the assembly of K r . 2-D problems are so cheap and 
Newton's method so robust that there is little to be gained over the employment of (83) and (85) as described 
here. But for 3-D problems with penalties, our experience indicates that the sensitivity of the iterative 
schemes and the dominant computational cost of assembling K r warrant the investigation of other iteration 
schemes. A likely candidate may be "modified Newton's method" in which K r is reassembled less often 
than once every step [21,22]. 



4.3. Interpolatory mesh refinement 

Another way to generate an initial guess either for the next load increment in (83) or to initiate (85) is to 
take a solution to (24) or (25) which is essentially fully converged on a given mesh and interpolate that solu- 
tion to a finer mesh. This is done by taking the values given by (12) at the new nodal locations to define 
nodal values of a new trial function in the form (12). Of course, the new trial function has more degrees of 
freedom, so in general it will not represent an equilibrium solution to (24) or (25). The hope is that it will be a 
good initial guess to such a solution. A strategy which has proved successful in our problems is to combine 
load increments (83) and any required refinements by Newton's method (85) on a crude mesh until A = 1. At 
the final load-level we refine the mesh in several steps — each generating an initial guess for the next which 
is refined by Newton's method. We have not performed a critical comparison of this strategy with others 
which are possible. The interpolatory refinements could be mixed with iterations of (83) or (85) in a variety 
of ways, and it would take careful study to determine the most efficient strategy. The reason we have used 
this particular strategy is that it provides a sequence of refinements at the final load-level which can give 
some idea of convergence of the model as the mesh is refined. 

Again we find that the 3-D model is much more sensitive to the mesh refinement scheme than the 2-D 
model, because of the penalty. The sensitivity is reflected in the fact that the interpolated initial guess is 
often outside of the region of convexity. As described in [10], this can be traced to the fact that volumetric 
integration points, £ r , of (25) and (29) on the refined mesh do not overlay the £ r on the cruder mesh. This 
leads to inaccurate computation of the pressure, h, via (41). A smoothing technique for h is described in [10], 
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which is based on the equivalence theorem of [8]. It is employed as a standard feature of the 3-D code and 
appears to work as long as the refinement is not too drastic. 

4.4. The in-core solvers 

The first priority in the choice of linear equation-solving procedures for the systems in (83) or (85) is that 
they take advantage of the symmetric, band-matrix structure of K r , described in (30). We store K T as d(co + 
1) vectors of length dN. Since it is hoped that the application of the procedures outlined earlier in this sec- 
tion will lead to positive-definite matrices, we employ elimination procedures with no interchanges 
[23,27-29]. There are several algorithms which can reduce the matrix to a compact triangular form which 
overlays the original matrix bands and require no additional storage space [23,27-29]. We employ two such 
algorithms. The first is the Cholesky method [23,27-29]; this will only work for positive-definite matrices. If 
any iterate of (83) or (85) departs from the convex region, the Cholesky method will fail, the iterations will 
cease, and our program will produce a suitable error message. 

The Cholesky method is used unless it is determined that there is no reasonable load-refinement strategy 
to avoid departure from the convex region. In that case we employ the symmetric column elimination pro- 
cedure described in [23]. This will usually provide the desired triangular decomposition, though, strictly 
speaking, the existence of the decomposition is not guaranteed without the use of interchanges [23,28]. 
Interchanges are not allowed if the band structure of the decomposition is to be preserved, and violation of 
band structure is incompatible with the goal of maximizing the size of K r which can be decomposed in-core. 
Even when the decomposition without interchanges exists for an indefinite K r , the decomposition process is 
potentially unstable in a numerical sense [28]. When this difficulty is combined with the tendency of 
Newton's method to become lost if convexity is lost, it means that the symmetric elimination must be used 
with caution and as a last resort, when convexity cannot be maintained. If the iteration procedures (83) and 
(85) survive the loss of convexity, and the iterates converge, eventually convexity must be regained [26]. 
When this happens, then there is no essential difference numerically or theoretically between the symmetric 
elimination and the Cholesky decomposition, as long as convexity is maintained. We find that in 3-D prob- 
lems, we often cannot avoid loss of convexity in the iteration scheme. On the other hand, if the load incre- 
ments, AA, are not too large and the mesh refinements not too drastic, then the loss of convexity can be 
survived and solutions obtained. Therefore the symmetric elimination provides an indispensible way of 
recovering from loss of convexity, which would be impossible with the Cholesky method. 

5. The stretching of rubber sheets 

To illustrate the type of problem the model described here is designed to solve, we give an example of the 
stretching of a rubber sheet This is similar to a problem described in [21], except that we have used the 
geometry of an A.S.T.M. standard tensile-test specimen for rubber sheet [30]. Our major concern was to 
investigate the application of penalty methods to impose the incompressibility constraint; therefore the 
emphasis in the results reported here is on the success or failure of the modelling process. However, we have 
been able to make some observations — largely qualitative — about the nature of the strain field in such 
specimens. 

5. 1 . The A.S.T.M. tensile-test specimen 

Figure 7 illustrates a top view of the A.S.T.M. tensile-test specimen described in [30] for the testing of the 
failure strength of rubber sheets. The dimensions are given in the lower portion of the figure. These speci- 
mens are referred to as "dumbbells" by the A.S.T.M., but this does not imply axisymmetry. Rather, the 
specimens are cut from a sheet which is typically one to several millimeters thick. Therefore, the plane-stress 
assumption described earlier seems appropriate. The 3-D model assumes a thickness of 2.5 mm. 

In the upper portion of figure 7, a fairly crude irregular mesh is depicted. Experimentation with various 
mesh subdivisions allows concentration of subdivision in areas where spatial variations of strain is the larg- 
est The large elements are placed in regions where the deformation is close to simple extension, which can 
be represented exactly, even by large elements of the type described in section 3.3. 
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FIGURE 7. The A.S.T.M. tensile test specimen, lower 
half showing fl with exact dimensions. Upper half 
shows a typical non-uniform finite element mesh giv- 
ing fl„. 
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5.2. Two- and three-dimensional models 

a. Symmetry 

As has been illustrated in figures 5 and 6, there are certain symmetry assumptions in the current imple- 
mentation of our model. Both 2- and 3-D versions assume symmetry of the deformation between the upper 
and lower portions of figure 7, as implied by figures 5 and 6. Also as implied by figures 5 and 6, symmetry of 
solutions between right and left halves of figure 7 is assumed, so that in the 2-D case, only the upper right 
quadrant is discretized. In the 3-D model symmetry is assumed along the midplane as implied by figure 6, so 
that only one octant of the body is actually discretized. 

b. Boundary conditions 

The boundary conditions along symmetry lines and planes are enforced as they are in linear analysis 
[22,23]. Displacements into, and certain shears along the symmetry lines and planes are assumed to be zero. 
The physical boundary conditions are a fairly simplified idealization of possible attachments to physical 
gripping devices (for example, see [31]). The conditions amount to a simple specification of end displace- 
ments to give various total extensions of the specimen. In 3-D, this roughly corresponds to gluing the whole 
end of the specimen to a relatively rigid surface. In the 2-D model, since the only end displacement available 
is on the midplane, the 2-D boundary condition is roughly like gluing the centerline of the end of the speci- 
men to a rigid surface. The rest of the end is free to slip on the support The 2- and 3-D physical boundary 
conditions are illustrated in figure 8. 

More realistic boundary conditions and interfacing with gripping devices can be envisaged, and could be 
incorporated in the model. In [32] methods for dealing with more realistic grips are described. However, for 
the purposes of testing penalty methods we did not feel that such considerations were crucially important 



BONDED SURFACE 




Figure 8. Boundary conditions for 3-D (upper)— bonding to a rigid support, and for 2-D flower)— bonding of 
centerline only. 

c. The use of penalties 

The penalty method of section 2.2 was employed to impose near incompressibility. The Mooney-Rivlin law 
(52) was used with Cl /c, = 0.2 and * = 50 c,. The actual value of c, serves only to scale the stresses, so c, can 
be used to give appropriate dimensional units. This value of the penalty leads to a typical compressibility 
error of about 2-4 percent For example, the most highly refined 3-D mesh of figure 7 has three elements 
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across the half-width and fifteen along the length and one through the thickness (3x15x1). For the element 
in this mesh with the highest strain component, III = 1.0344 instead of III = 1. 

A penalty is also used to enforce the displacement boundary condition on the end of the specimen [14]. If 
Va denotes the surface of the specimen which is bonded to the support, then the following term is appended 
to the variational principle (6): 

z d \ Td (u i -u*)(ii i -u,*)ds (sum on i) (86) 

where z d is the boundary penalty and «? is the specified displacement vector. In practice it has been found 
that z d can be taken quite large. In the current examples we have taken z d = 10 4 . The boundary integral in 
(86) is easy to evaluate [22,23], Employment of this boundary penalty means that all essential b.c. are 
homogenous [22,23]. 



5.3. Code verification procedures 

Because exact solutions to the problem we describe here are not available, special attention must be paid 
to insuring that the computer code is functioning correctly. Of course, it is not really possible to give an 
ironclad guarantee of the correctness of the code, but we have used three tests which assure us, beyond 
reasonable doubt, that the results we present here are artifacts of the model and not coding blunders. 

a. The small-strain limit 

As described in section 3.4, our code has the capability to produce linearly elastic solutions with the same 
geometry as for the nonlinear results. A simple asymptotic analysis of the Mooney-Rivlin law and the varia- 
tional principle (6) using that law shows that when squares of the components of G y — d, 7 are small enough to 
be neglected, linear elasticity and nonlinear elasticity should give nearly the same results. We found that 
our nonlinear code reproduced the 2- and 3-D linear displacements to 4 digits when the total extension of the 
specimen was on the order of 0.5 percent. The differences between the two solutions were typically in the 
second digit for extensions in the range of 1-2 percent and progressively more noticeable for higher exten- 
sions. These results are in total accordance with what is to be expected from the asymptotic analysis. 

b. Quadratic convergence of Newton's method 

The fact that the small-strain limit is obtained from the nonlinear code does not guarantee that the code is 
working correctly. If an error is made in the coding of energy-density terms of higher order in displacement- 
gradients, this might not affect the small-strain limit In the course of debugging this code and other 
nonlinear codes, it has been our experience that the quadratic convergence of Newton's method described 
in section 4.2 can only be obtained if K r (u) is derived and coded correctly. What the obtaining of quadratic 
convergence means is that the derivatives computed according to section 2.3 are the correct second 
derivatives (29) of the first derivatives (24) of the variational principle (6). It does not completely assure that 
the first derivatives are derived correctly. The fact that, when solutions are obtained, K T (u) is invariably 
positive definite suggests that the first derivatives are also correct. 

Table 1 shows the obtaining of quadratic convergence in a 2-D mesh (pt. A.) and a 3-D mesh (pt. B.). Part 
C. of table 1 shows an example of a sequence of iterations for a 3-D mesh becoming "lost." All values in the 
table are the square roots of the sum of the squares of the components of P t . The 2-D results are from a 
7 X 26 element mesh discretizing the first quadrant of the specimen, extended to a total extension of 200 
percent Ten load increments (83) of 20 percent were performed on a regular 3x5 mesh with no intermedi- 
ate refinement Then the 3x5 fully extended displacements were refined by Newton's method and interpo- 
lated to a 4x 7 regular mesh. The same procedure was followed to proceed to a 5 X 10 mesh, a 5 X 12 mesh, 
then a 6x 15 mesh, and finally these results were interpolated to the irregular 7 X 26 mesh. It is the refine- 
ment by Newton's method of these last interpolated results which is given in table 1A. 
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TABLE 1. Quadratic Convergence of Newton's method as reflected in the decline of the rms residual of the equilibrium equations 

A and B. C shows iterations becoming lost. 

A. B. C. 

iter, no. 7 X 26 2-D 3 X 15 X 1 3-D 3x8x1 3-D 

1 0.422 X 10 1 0.551 X 10 3 0.900 X 10 3 

2 0.588 X 10° 0.271 X 10 2 0.794 X 10 2 

3 0.262 X 10" 1 0.665 x 10 1 0.529 X 10 s 

4 0.652 X 10" 4 0.883 X 10° 0.370 X 10" 

5 0.627 X 10"' 0.108 X 10° 0.130 X 10« 

6 0.445 X 10" 14 0.585 X 10" 3 0.176 x 10" 



Table 1 B is obtained from the performance of 10 load increments of 10 percent on a 3 X 8 X 1 element 3-D 
mesh, to give a total extension of 100 percent Newton's method was used to refine the results after every 
two load increments. Interpolatory refinement with Newton iterations was used to go to a 3 X 10 X 1 and then 
a 3 X 12 X 1 regular mesh. Finally, the results were interpolated to the irregular 3x15x1 mesh whose top 
view is given in figure 7. It is the refinement by Newton's method of these last interpolated results which is 
given in table IB. The LDL r column solver is required for the procedures involved in table IB, since con- 
vexity is lost after the first pair of load increments. It is always regained by the Newton iteration, but is lost 
at the next load increments, and just after each interpolatory refinement. 

Table 1C shows what happens when an attempt is made to carry out the load increments described for 
table IB on a 3x5x1 mesh instead of a 3x8x1 mesh. The loss of convexity is so severe that when an 
attempt is made to interpolate the displacements to a 3x8x1 mesh and refine by Newton's method the 
results of table 1C are obtained. The ability to obtain quadratic convergence with a careful choice of load 
increments leads us to believe that the code is correct, but that remaining close enough to the convex region 
is a delicate matter in the 3-D model. 

c. Convergence with mesh refinement 

Figure 9 shows values of the extensional strain J 12 and shear strain J 2l for a specimen extended to a total 
extension of 100 percent. The results are taken from a 2-D calculation with 26 elements along the length, 
with smaller elements concentrated at the top of the neck. The smallest elements are on the order of a 
quarter the size of the smallest elements in figure 7 and the larger elements of the refined 7 X 26 mesh are in 
the simple extension portion of the neck and are half the length of the corresponding elements in figure 7. 
Figure 9 plots values of J 22 and J 2l at a position indicated by the arrow on the cruder 3x 15 mesh (both 
halves shown) of figure 10. This is where the strain-gradients are largest. With the same lengthwise subdivi- 
sion of the specimen into 26 pieces, results are compared for meshes which have 3, 4, 6 and 7 elements 
across the neck. The 7x26 mesh is obtained from the 6x26 by subdividing the outer element in two. Also 
plotted on figure 9 are the results of the 3x15x1 3-D mesh, whose top view is given in figure 7. 

The results show that to a large extent, the results from all the meshes overlay the same curve, with the 
largest deviations coming — as expected-from the 3x 15 X 1 and 3x26 meshes. Therefore we are led to con- 
clude that convergence is taking place and that, qualitatively at least, even the cruder meshes give good 
results to graphical accuracy. A similar conclusion can be drawn from a comparison of tables 2 and 3 in 
which values of the pressure H of (70) for a 7x26 and 3x 15 mesh are compared — likewise for a total exten- 
sion of 100 percent 

We believe that the circumstantial evidence gathered from the obtaining of the small-strain limit quad- 
ratic convergence of Newton's method, and convergence with mesh refinement make it highly unlikely that 
our numerical results are an artifact of coding blunders. We have run exactly parallel two- and three- 
dimensional versions of the 3x 15 mesh of figure 7 (3x 15 X 1 in the 3-D case). While we have not run as 
many 3-D convergence tests as 2-D convergence tests, we find that the agreement in displacements and 
strain is in general as good as the agreement illustrated in figure 9. We are convinced that convergence of 
displacements with mesh refinement is taking place in the 3-D model. What is happening to the pressures in 
table 4 is another — and very interesting — matter, which will be discussed below. 
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FIGURE 9. Extensional strain (upper) and shear-strain (lower) sampled at element centroids in stress-concentration region for various 
meshes. Values given on cross-section of neck. 




FlGURElO. Mesh of figure 7 deformed to total extension of 100 percent Arrow marks stress- 
concentration. 



Table 2. Piecewise constant pressures along a cross-section of the neck in the stress-concentration region. Values taken from a refined 
2-D mesh (7 * 26). Positions taken at center point of individual elements. 



Position 
Pressure 



0.25 
0.537 



0.75 
0.533 



1.25 
0.523 



1.75 
0.508 



2.25 
0.480 



2.625 
0.452 



2.875 
0.438 



TABLE 3. Similar values to table 2, taken from a crude 2-D mesh. 



Position 
Pressure 



0.50 
0.534 



1.50 
0.514 



2.50 
0.498 



Average value 
0.515 



TABLE 4. Values of the pressure taken from a 3-D mesh along cross-section near the 
stress-concentration region. 



Position 
Pressure 



0.50 
0.621 



1.50 
0.454 



2.50 
0.563 



Average value 
0.546 
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5.4. Discussion of results 



a. The plane-stress results 



Our experience shows that the use of Newton's method with incremental loading described in section 4., 
when combined with the isoparametric plane-stress model described in this paper, produces a robust and 
computationally inexpensive algorithm. However, we must point out that the geometric assumptions and 
boundary conditions used in the current study advise a cautious interpretation of our results. We do think 
that we can draw some broad conclusions about the behavior of the A.S.T.M. speciments based on the quali- 
tative picture which emerges from our 2-D results. 

The picture which emerges is that the A.S.T.M. specimen geometry does what it is intended to do: namely 
concentrate the extension of the specimen in the neck region and away from the grips. The geometry does 
produce a region of concentrated deformation at the position of the arrow in figure 10, but our results indi- 
cate that the concentration is not severe. The data we cite to support these points, some of which are plotted 
in figure 9, are taken from our 7 x 26 mesh — not the cruder mesh of figure 10. 

To illustrate the concentration of extension in the neck, we look at the value of the extension ratio, J 22 , in 
the element nearest the center of the neck as compared to the total extension ratio of the specimen (based on 
grip separation). For 10 percent total extension of the specimen, the neck value of J 22 is 7.2 percent higher 
than the total extension ratio of 1.10. For 100 percent total extension, the neck has aj 22 51 percent higher 
than the total extension ratio 2.00, and for total extension of 200 percent, J 22 is 67 percent higher on the 
neck than the total extension ratio of 3.00. 

Turning to the region of concentrated deformation, we find that there is a maximum of J 22 at the point 
indicated by the arrow in figure 10. For 10 percent total extension, the maximum extension ratio is 5.5 
percent higher than the extension ratio in the neck. For 100 percent total extension, the maximum extension 
ratio is 10 percent higher than in the neck, and for 200 percent total extension, the maximum is 20 percent 
higher than in the neck. 

The extension ratio in the region of concentrated deformation seems to approach a value of 20 percent 
above the neck value in the limit as total extension is increased. The extension ratio J 12 is not a rotationally 
invariant quantity, and one might raise the question as to whether the concentration of deformation grows 
unboundedly with total extension ratio by increasing shear-strain concentration. We look at the energy- 
density, which from (50), (51) and (64) can be seen to be W(l, II) = c t (1-3) + c 2 (II— 3), if one requires that 
the energy-density of the undeformed state is zero. Again we compare neck values of W with the maximum 
value. For 10 percent total extension W is 92 percent higher at its maximum than in the neck; for 100 
percent total extension, it is 57 percent higher; and for 200 percent total extension, it is only 53 percent 
higher. So the deformation-energy concentration actually decreases in intensity with increasing extension 
ratio. Thus the shear-strain intensity cannot be growing unboundedly. In fact, if we look at the quantity 
J 2 i/(J 2 2~1), we find that for 10 percent total extension, the maximum value is 27 percent above the value 
in the neck. For 100 percent total extension, the intensity is 21 percent over the neck value, and for 200 
percent total extension, the intensity is 19 percent over the neck value. So in fact shear-deformation appears 
to become relatively less important as the total extension increases. 

The A.S.T.M. standard [30] specifies that the total extension of the neck be measured visually, using 
marks on the neck. If the results we present here represent an adequate qualitative picture of the specimen, 
we would expect the specimen to fail in the neck, just below the wide portion. The A.S.T.M. standard 
requires that neck extension be recorded at the moment of failure. Our results indicate that there is a built- 
in safety factor of about 20 percent, in that the actual extension ratio near the point of faliure could have 
been about 20 percent higher than recorded. Shear deformation near the point of such a failure would not 
appear to be a major factor. 

b. The three-dimensional results 

For reasons discussed in detail earlier, the 3-D incremental loading/Newton algorithm combined with the 
isoparametric finite element/penalty model described here cannot be classified as a robust or computa- 
tionally inexpensive algorithm. One would expect 3-D nonlinear finite element models to be dramatically 
more expensive than 2-D models in general But one would also hope that the nonlinear iteration scheme 
could be made more reliable, so that situations like the one illustrated in table 1C do not lead to the waste 
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of expensive iterations. But there is a more basic theoretical question raised by the 3-D results, which 
should be answered before much effort is spent on the optimization of iteration schemes. This is a question 
of the stability and accuracy of the penalty or equivalent mixed formulation [8, 13-16]. 

Table 5 shows values of the pressure in a cross-section of elements near the stress-concentration region of 
a 3 X 10 X 1 mesh. The elements are rather large, and the centroids from which the pressures are taken 
are further down the neck towards the middle of the specimen than the sample points of tables 2 and 3, 
therefore the pressure is lower. But note that there is a smooth variation of pressure. Ths contrasts sharply 
with the oscillation of pressure across the cross-section in the stress-concentration region shown in table 4. 
These results are taken from the irregular 3 X 15 X 1 mesh whose top view is given in figure 7. We have 
already seen that J 22 and J 2i vary smoothly and are reasonably accurate in these elements from figure 9. 
In general, a detailed comparison of the 3 X 15 X 1 results shows that u t and/ y compare well with those 
from the 2-D 7 X 26 mesh — particularly in comparison to the contrast of the results in tables 2 and 4. What 
we observe there is that the 3-D pressures alone are affected by a "checkerboard" pattern of high and low 
values in the stress-concentration region. This is illustrated in figures 11 and 12. The pattern of high-low- 
high on one cross-section is countered by a reverse pattern on adjacent cross-sections. The average pressure 
value along a checkerboarded cross-section seems to agree much better with average values taken from the 
2-D model, as is illustrated in figure 12 and tables 3 and 4. This leads us to speculate that averaging 
techniques may be able to remove the checkerboard. 



TABLE 5. Values of the pressure taken from a cross-section near the stress 
concentration for a cruder 3-D mesh. 



Position 
Pressure 



0.50 
0.446 



1.50 
0.443 



2.50 
0.421 




no. 4 5 



FicUREll. Top view of a 3-D mesh— qualitative picture of 
checkerboard pressures. 
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FiGUREl2. Piecewise-constant pressures on a cross-section of the neck: (a) Accurate pressures; (b) 
checkerboard pressures; (c) checkerboard variation on adjacent cross-section to(b). 



In view of the success obtained using this trilinear isoparametric element in a variety of penalty or mixed 
formulations [2, 7, 8, 11], we are led to conclude that the checkerboard pressure problem is a relatively rare 
occurrence. But, particularly since this phenomenon occurs when the mesh is refined, we are led to doubt 
that our element can satisfy the basic stability requirement of [13-16]. This matter is under investigation 
by the first author and others [14, 15, 34]. Techniques to overcome these difficulties using modified 
penalty procedures are also being investigated [18, 19, 33, 34]. Though the appearance of checkerboard 
pressure keeps us from drawing any physical conclusions based on our 3-D results, we are convinced that 
it will prove extremely useful to future theoretical investigations to be able to reliably reproduce this 
checkerboard phenomenon in such a basically simple problem. 



6. Conclusions 

The purpose in the development of the finite element model described here is to generate test problems 
for research on penalty methods in finite elasticity. As such it has proved very successful. In spite of the 
simplicity of the model, it has been able to isolate two important areas which will require further investi- 
gation if finite element penalty methods are to become productive tools in the analysis of incompressible 
materials under finite strain. The first area is that of the sensitivity of iteration schemes to small volume 
changes. Progress has been made in this direction, but still more is required [9, 10]. Perhaps more important 
is the question of element stability raised by the appearance of checkerboard pressure modes. The 
stability condition is known [13-17], but few elements seem to satisfy it— whether in penalty or equivalent 
mixed formulations [8-10]. On the other hand, even for elements which apparently do not satisfy it, the 
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checkerboard modes may be removable and often do not occur. More precise mathematical predictions 
need to be provided for these phenomena. It would be very useful to be able to aid such theoretical 
investigations with a code with more general element capabilities than the current one. Clearly the model 
admits to a variety of element types. Beyond that, new developments indicate that there are other penalty 
formulations which may be more stable than the reduced/selective integration techniques incorporated 
in the current model [19, 33, 34], Therefore the model should be enhanced to include these possibilities. 

Finally it is tempting to think of ways in which to make the code itself more sophisticated. The mesh 
generator could be based on the unions of several reference bricks, making the variety of bodies describable 
much wider. Load increment/ refinement strategies could be made adaptive, in order to maintain convexity 
as much as possible without human intervention. More physically realistic boundary conditions and non- 
conservative loads could be incorporated. Many other constitutive laws fit within the framework of the 
model and would require only minor coding to implement However, these enhancements belong further 
down the road, after some of the challenging questions raised by the current model are resolved. 



The finite element model and the computer code described here were developed while the first author 
was a Postdoctoral Research Associate at NBS. He would like to express his appreciation to that organiza- 
tion for the opportunity to carry on this research, which continues to be a fruitful avenue of investigation. 
He would like to express his gratitude to his advisors, E.A. Kearsley and J.T. Fong and thank Applied 
Mathematics Division Chief, B. Colvin, for his support and encouragement 
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